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A BSTRAC T 

Oiten in the short-range limited-area numerical weather prediction 
(NWP) of extfdtropicel weather systems the effects of planetary boundary 
layer (PBL) processes are considered secondarily important. However, it 
may not be the case for the regional NWP of mesoscale convective systems 
over the arid and semi -arid highlands of the southwestern and south-central 
United States in late spring end summer. Over these dry regions, the PBL 
con grow quite high up into the lower middle troposphere (600 mb) due to 
very effective solar heating and hence a vigorous air-land thermal 
interaction con occur The interaction representing a major heat source for 
regional dynamical systems can not be ignored. 

The present study focuses on the development of an one-dimensional 
PBL-rodiation model. The model PBL consists of a constant-flux surface 
layer superposed with a well-mixed (Ekman) layer. The vertical eddy mixing 
coefficients for heat and momentum in the surface layer ere determined 
according to the surface similarity theory, while their vertical profiles in 
the Ekman layer are specified with a cubic polynomial. Prognostic equations 
ore used for predicting the height of the nonneutrel PBL. The atmospheric 
radiation is parameterized to define the surface heat source/sink for the 
growth and decay of the PBL. A series of real-data numerical experiments 
has been carried out to obtain a physical understanding how the model 
performs under various atmospheric and surface conditions. 

This one-dimensional model will eventually be incorporated into a 
mesoscale prediction system. The ultimate goal of this research is to 
improve the NWP of mesoscale convective storms overland. 
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I. INTRODUCTION 


Often in the short-range limited-area numerical weather prediction 
(NWP) of extratropical weather systems the effects of planetary boundary 
layer (PBL) processes are considered secondarily important. This is 
probably true over land in winter or over the oceons in summer when the PBL 
is generally very shallow. However, it may not be the cose for the regional 
NWP of mesoscale convective systems over the arid and semi-arid highlands 
of the southwestern and south-central United Stotes. Over these arid 
regions, in late spring and summer the PBL can grow quite high up into the 
lower middle troposphere (600 mb) due to very effective solar heating. This 
can cause vigorous air-land thermal interaction. The interaction 
representing a very significant heat source/sink for regional dynamical 
systems can not be ignored. 

The present study will focus on the development of an one-dimensional 
model coupling the PBL processes and atmospheric radiative transfer. The 
model PBL consists of a constant-flux surface layer superposed with a 
well-mixed (Ekman) layer. The vertical eddy (turbulence) mixing 
coefficients for heat and momentum in the surface layer will be determined 
according to the surface similarity theory, while their vertical profiles in 
the Ekman layer are specified with a cubic polynomial. Prognostic equations 
are formulated for predicting the height of nonneutral PBL. The atmospheric 
radiation including solar and infrared components is parameterized to 
define the surface heat source/sink for the growth ond decay of the PBL. 

This one-dimensional model will eventually be incorporated into the 
limited-area mesoscale prediction system (LAMPS) (Perkey,1976; Chang et 
ol. , 1 98 1 ) for regional NWP over land. 

The model levels range from 0 to 16 km with relatively higher 
resolution in the PBL. A discussion of PBL physics ond techniques for 
solving the atmospheric radiative transfer in this study is included in 
Sections 3 ond 4, respectively. The results of three numerical experiments 
conducted on the PBL-rodiotion model ore described in Section 5. 


2. OBJECTIVES 

The ultimate goal of this research is to improve the NWP of mesoscale 
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convective storms over arid land. Before incorporating the PBL-radiation 
model into LAMPS, we like to carry out a series of real-date numerical 
experiments to have a physical understanding how the one-dimensional 
model performs under various atmospheric conditions. The specific 
objectives include: 

1) To assess the impact of clouds on the surface heat sources and on 
the PBL development. 

2) To determine the sensitivity of the PBL height to its initial value. 

3) To understand the impact of surface conditions such as roughness 
length, albedo, the Bowen ratio on the structure of the PBL. 

4) To understand the role of the static stability in the growth and 
decay of the PBL. 


3. PLANETARY BOUNDARY LAVER 

The PBL is a layer of atmosphere on the order of 1 km in depth above the 
earth’s surface. This thin layer of air is characterized by small-scale 
turbulence of spatial dimension no greater than 1 km. In the free 
atmosphere above the PBL, the turbulent motions are considerably weaker. 
The PBL processes represent a consequence of interaction between the 
lowest layer of air and the underlying surface. The interaction can 
significant impact on the dynamics of the upper air flows. 

The influences of the small-scale eddies on large-scale (model 
resolable scale) atmospheric circulations moy be included in the model 
equations as described below. These equations ore formulated on the 
(x,y,z,t) coordinate system, in which x is along the lotitudinol (west to 
east), y is along the longitudinal (south to south) circles on the earth’s 
surface, while z is perpendicular to the x-y plane. 

3.1 The basic equations 

Equations of horizontal motion. 

du/dt = fv - p -1 dp fbx - p'Hdp^'uVdx + dp<u'v’>/dy ♦ dp<u‘w’>/dz] (1) 

dv/dt =-fu - p' ! dp/dy - p'Hdp^’vVdx + dp<v'v’>/dy + dp<v M w’>/dz) (2) 
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Hydrostatic equation 

dp/dz * -pg (3) 

Equation of state for air: 

p = pR0(p/1OOO)° 28571 (4) 


Mass conservation: 

dp/dt = -p(du/dx + dv/dy + dw/dz ) (5) 

Energy conservation: 

d0/dt = H + p' , [d<-pu'0'>/dx + dp<-v‘6'>/dy + dp<-w'0'>/dz] (6) 

Moisture conservation: 

dq/dt = Q + p" , [d<-pu'q>/dx + dp<-v'q’>/dy + dp<-w’q’>/dz] (7) 

where d /dt = d /dt + ud /dx ♦ vd /dy + wd /dz. There ore seven basic 
large-scale variables in the above system: three velocity components in the 
x, y, and z directions, respectively, u ,v, and w; pressure p; air density p; 

potential temperature 0, and specific humidity q. Other symbols include the 
Coriolis parameter f, heat source H, moisture source Q, and the gas constant 
for dry air R. 

The prime ' denotes the small-scale turbulent fluctuations, while the 
bracket < > denotes a time average over a period much longer than the time 
scale of the turbulent fluctuations. In the u- and v-equations, the first 
term represents the Coriolis force due to the earth's rotation, and the 
second term represents the pressure gradient force. The six covariances 
(e g., -p<u u >) are referred to as eddy stresses representing eddy or 
turbulent fluxes of momentum. Similarly, there are eddy flux of heat in the 
0 equation and eddy flux of moisture in the q equation. 

3.2 Closure problem 
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Considering the basic variables only, we note that the above set of 
equations forms a closed system. Given initial conditions, the future states 
of the seven basic unknowns can be determined in principle by numerical 
integration. However, with the eddy flux terms representing additional 
dependent variables additional equations are required to close the system. 
This is so-called closure problem in dealing with atmospheric turbulence. 

In the free atmosphere, for the atmospheric circulations of horizontal 
scale larger than 100 km and time scale on the order of one day or longer the 
contribution to the local change (d /dt) due to the eddy stresses as well as 
the eddy fluxes of heat and moisture is insiginificant in comparison with 
those due to large-scale processes (e g., pressure gradient force, advection, 
etc ). In NWP models, we may neglect the eddy flux terms in the above set 
of equations. This is known as the zero-order closure. 

In the PBL, however, the eddy flux terms may have the same order of 
magnitude as the other terms and can not be neglected under most 
situations. To close the system without neglecting PBL processes in a NWP 
model, we may consider the following two closure approaches. 

a) First-order closure 

A well-known and commonly used technique in PBL meteorology for the 
first-order closure is the so-called K-theory. The theory simply states that 
the turbulent flux of any physical quanties such as momentum and heat may 
be treated analogously to molecular diffusion. For example, for the vertical 
eddy flux of heat, we may write 

;0‘w‘> = Kd0/dz 

which indicates the vertical mixing of heat as a result of turbulence is 
parameterized in terms of the large-scale variables. Similar expressions 
can be formulated for the eddy mixing of momentum. A more detailed 
discussion of the K-theory and the procedures to implement the theory to 
PBL modeling in this research is presented in Section 3.3. 

b) Higher-order closure 

The concept consists of deriving prognostic equations for the eddy flux 
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terms from the bosic equations. The derived equotions contoin higher-order 
correlations. For example, in the case of second-order closure triple 
correlations such os <u‘u'v‘> will appear in the prediction equation for <u'u >. 
Thus, the closure problem remains and closure assumptions have to be made 
at a higher level to complete the system. Moreover, the introduction of the 
equations for the eddy components in NWP models will drastically increase 
the demand for computer time. Also, the current observations are 
insufficient to provide an accurate initial condition for the eddy fluxes. To 
date, higher-order closure approaches are not yet desirable for NWP models. 

3.3 First-order closure PBL model 

As shown in the basic equations, the eddy mixing of momentum, heat, 
and moisture includes both horizontal and vertical components. Due to the 
horizontally stratified nature of the atmosphere, the gradient of u, v, 8, and 
q is much greater in the vertical than in the horizontal direction. The eddy 
mixing in the horizontal is negligible comparing to its counterpart in the 
vertical. Hence, with the use of the K-theory the central issues in PBL 
research involve the determination of the vertical profiles of K for 
momentum, heat, and moisture as well as the depth of the PBL within which 
the profiles are applied 

In formulating equations for the K profiles, a two-layer model is 
developed. The PBL is devided into a constant-flux surface layer superposed 
with a well-mixed Ekman layer. 

a) Surface layer 

Tiiis is a very thin turbulent layer within a depth less than 50 m above 
the surface. In this layer, the turbulent fluxes are nearly constant with 
height and the vertical wind profile depends on the static stability and 
surface roughness. There are mechanical turbulence driven by the vertical 
wind shear as well as thermal turbulence driven by buoyancy (dry 
convection). The former decreases rapidly with height, while the latter 
varies slowly with height. 

For a surface layer of neutral or near-neutral stability (<8'w'> * 0 or 
mechanical turbulence only, the vertical wind profile is very close to the 
well-known logarithmic form 
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and the corresponding K, eddy viscosity for both u and v, is 


K - kzu* (9) 

where k = 0 35 the von Karmon constant, u* the friction velocity, and z 0 the 

roughness length. Over a given surface, z 0 is a known constant, while u* can 

be evaluated from the observed wind profile. Thus, K is determined. 

Note that K increases upward suggesting the eddy size is proportional to the 
height above the surface. 

For a nonneutral surface layer, there will be vertical heat exchanges 
between air and the underlying surface. We will have predominant 
mechanical turbulence if the layer is stable 08/dz > 0), or both mechanical 
and thermal turbulence if the layer is unstable O0/dz < 0). According to the 
similarity theory summarized in the book by Haltiner and Williams (1980), 
the vertical gradients of the wind and potential temperature can be 
expressed as follows, respectively. 

du/di = [(u#/k)/z] 4> m (z/L) (10) 

dQ/dz = [(0*/k)/z] cf?j-j(z/L) (11) 

L = 8 o u* 2 /(kg0*) (12) 

where e* analogous to u* representing a temperature scale related to the 
vertical eddy heat flux, and 0 O the potential temperature at the surface. 

The length scale L, called the Mon in-Obukhov length derived based upon 
dimension analysis, measures the relative role between mechanical and 
thermal forcing in generating the turbulence. L is negative in unstable air 
with upward heat flux, while positive in stable air with downward heat flux. 
Also, physically the absolute value of L may be interpreted as the height 
above which the mechanical turbulence becomes insignificant compared to 
the thermal turbulence. Thus, as air becomes more unstable, the absolute 
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value of L becomes smaller. <J> m and <}> h are two nondimensional universal 
functions determined empirically from observations. The corresponding 


eddy viscosity K m and eddy conductivity K h are 

K m = ku*z/<t> m (13) 

K h = ku«z/<|> h (14) 

For unstable air 

<J) m « (1 - 15z/L]- ,/4 and 4> h = 0.74[1 - 9z/L]' |/2 (15) 

For stable air 

<J) m * [1 ♦ 4.7Z/L] and - 0.74 + 4.7z/L (16) 


For moisture flux, K q known as eddy diffusivity is assumed to be the same 
as K h . 

b) The Ekman layer 

The vertical extent of the PBL is variable. Its depth may range from 100 
m (e g., the nocturnal PBL) in stable air to over I km (eg., the daytime PBL 
over arid land) in unstable air. Out of the total depth the surface layer 
occupies less than 10% of the lowest part, and the rest is the well-mixed 
Ekman layer. Because of its elevation, it is much more difficult to conduct 
observational study in the Ekman layer than in the surface layer. As a 
result, there is lack of a simple and elegant theory, such as the surface 
similar theory, which may be used to define the K profiles in this upper part 
of the PBL. Here, a cubic polynomial similar to the O'Brien’s formulation 
( 1970 ) except for the constraint of zero slope at the top of the PBL is 
adopted for the the K profiles in the Ekaman layer. 

The cubic polynomial has the following form: 
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( 17 ) 


K(z) - K(z s ) [l-9.75z + 8.75z 3 ] ♦ dK(z s )/dz [z e -z s ] [z-3.5z 2 +2.5z 3 ] ♦ 

K(z m ) [125z 2 (1-z)/12] ♦ K(z e ) [z 2 (5z-2)/3] 

where K(z s ), K(z e ), and K(z m ) are K evaluated at the top of the surface layer 
(z s ), at the top of the Ekman layer (z e ), and at the middle of the PBL (z m ), 
separately. K(z 5 ) and dK(z s )/dz are obtained from the surface layer physics 
discussed earlier, and K(z e ) has to be specified. Also, K has the maximum 
value at z m . In the current PBL model z s is set at 25 m and K(z e ) is set 
equal to I m 2 s''. 

c) Height of the PBL 

To apply the K profiles, we must know the depth of the PBL. Based on 
theoretical and observational studies, for the neutral PBL, Panofsky and 
Dutton (1984) suggested 

z e = 0.2u#/f (18) 

For the unstable PBL, the Deardorff's (1974) prognostic equation is used. 
dz e /dt = 1.8(w# 3 ♦ l.lu* 3 - 3.3u# 2 fz e )/(z e 2 o ♦ 9w« 2 + 7.2u« 2 ) (19) 

where w* = IgS^Zg/eQ)] 1 /2 , the vertical velocity scale; o = (g/0 o )d8/dz, the 

static stability at the surface. For the stable PBL, the following equation 
suggested by Nieuwstadt and Tennekes (1981) is used. 

az e /at = -udz e /dx -vdz e /dy + (z'-z e )/t s (20) 

where z‘= 0.4(u«L/f) |/2 , the equilibrium height; t^O.ysiQtegHQKdQQ/dt)' 1 , 
the time scale. 

d) Extremely unstable PBL 

Recently, Yoh (1989) suggested that under extremely unstable conditions 
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some modifications of the surface similarity theory were necessary. In 
very unstable air, -L approaches the order of meters and free convection 
responsible for generating large local eddies plays a key role in the vertical 
eddy mixing. It is inappropriate to use u#, which is strongly influenced by 

the mean wind, as a velocity scale to determine the eddy fluxes. Also, 
because <f> m (z/L) and fy^z/L) are accurate for a relatively small lz/Ll (< 5), 

both functions may only be applied at the level of a few meters in height. 

Yoh subdivided the surface layer into two layers: the lowest layer of 
1.22 m upward from the suface, in which the modified similarity theory will 
be applied, and 1.22-25 m layer, In which the mixed layer scaling will be 
applied. The major difference between the modified and the original 
similarity theory discussed earlier is the incorporation of w#, a free 

convection seal ing velocity, in computing the modified friction velocity, u + 

and Monin-Obuhkov length L+. The mixing length scaling technique is 

designed to relate u, v, 6, and q at 1.22 m to those at 25 m. This information 
provides an input for the modified surface similarity. The relationship 
between u + and L + is identical to that between u# and L (see Equ. 12). After 

L + is determined, Equations (13M16), with L replaced by L + , are used to 

obtain Ks at 25 m. Finally, the K profiles in the Ekman layer under very 
unstable conditions can be determined from (17). 

e) Surface heat source 

From the above discussion. It is clear that under the nonneutral 
conditions surface heat fluxes exert great influences on the K profiles in 
the PBL The intensity of these heat fluxes depends on to a large extent the 
temperature difference between air and the surface. To solve numerically 
the PBL equations, we often need tfte surface temperature as one of lower 
boundary conditions. 

Over the oceans, water responds slowly to atmospheric radiation and the 
temporal variation of water temperature Is small. A constant sea surface 
temperature (55T) field at its initial observed values Is adequate for model 
PBL in a short-range ( < 48 h) NWP. Over land, due to its quick response to 
solar heating the temporal variation of surface temperature is rather large. 



-v diurnal character is expected in the PBL height as well as the K profiles. 
Thus, a mechanism for determining the surface temperature and the 
associated heat flux is needed for properly resolving the essential features 
of the PBL over land This is particularly true over arid land. 

In this study, the land surface temperature will be calculated based on 
a heat budget consideration. The surface heat budget involves infrared and 
solar irradiance, and sensible and latent heat flux. Also, the surface heat 
budget is well known to be a strong function of cloud cover and moisture 
content in air In the next section, a brief description of the treatment of 
radiation processes as well os the procedures to obtain the surface 
temperature is presented. 


4 RADIATIVE ENERGY AND HEAT BUDGET AT THE SURFACE 

In modeling radiative transfers the usual constraints of a numerical 
model, i.e., not to be too time-consuming, must be considered. In this study 
only water vapor, which is one of the prognostic variables in the model 
equations, is regarded as a radiatively active gas in the model atmosphere. 
Ozone and carbon dioxide ore not considered in radiation computation 
because the model top is set at the lower stratosphere. In the troposphere, 
for a short-range NWP carbon dioxide contribution to the atmospheric 
radiation processes is negligible comparing with that due to water vapor. 

To solve the transfer equations for infrared radiation a simplified 
method based on the broad-band emissivity technique (Staley and Jurica, 
1970) is used. The technique has been successfully applied in many NWP 
studies (Danard, 1969, Chang, 1980). With respect to infrared radiation, 
clouds and the earth s surface are treated as a blackbody radiator. For 
modeling solar radiation, on empirical technique designed for the UCLA 
general circulation model (Haltiner and Williams, 1980) is used. Some 
modifications have been made to accommodate multi-layer clouds in the 
current model atmosphere The reflection of insolotion at the cloud top and 
the earth’s surface, and the absorption of insolotion by clouds are 
considered. The intensity of solar radiation varies with the zenith angle 
which is computed os a function of location and time. 

The sum of the insolation and the net infrared (downward minus 
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upward) radiation at the surface is taken as the total available surface heat 
source in the day or sink at night. The heat source/sink is partitioned into 
sensible heat flux, H s , and latent heat flux, Hj, according to the Bowen ratio 

(Sellers, 1967). The surface temperature is then related to sensible heat 
flux by a bulk aerodynamic formula of the form 

Hg = pCpC d [/(U 2 + V 2 )] (0 S - 0g) 

where c p specific heat capacity of air at constant pressure, 8 S the surface 
potential temperatures, and 8 a the potential temperature at 25 m. The drag 
coefficie.nce, Cp, is defined as 

Cp = k 2 [ln(z/z t ) - 4> h (z/L)]~’[ln(z/z 0 ) - <J? m (z/L)] -1 

where z t the roughness length for temperature. 


5. RESULTS 

Disregarding the large-scale processes in the basic equations shown 
earlier and with the use of K-theory for the small-scale eddies, we reduce 
the u, v, 0, and q equations to the following forms. 

du/dt = - p'^pK^u/dzydz 
dv/dt * - p _l d(pK m dv/dz)/dZ 
ae/at = - p-'a(pK h ae/az)/az 
dq/dt = - p-'a(pK h aq/dz)/dz 

These equations in conjunction with the hydrostatic equation and the 
equation of state for air constitute a closed system for the one-dimensional 
PBL-radiation model 
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The model hos15 levels in the vertical (z in Table 1). The PBL structure 
is resolved in terms of the lowest 6 levels. The time step used for the 
model integration was 1 minute and the radiation computation was updated 
every 30 minutes. The input sounding (Table 1) for the experiments was 
taken in Kansas at 1200 GMT on 1 1 June 1985 during the PRE-STORM. Other 
PBL parameters required for integrating the equations include the roughness 
length = 20 cm, the Bowen ratio = 0.05, and the surface albedo = 0.9. Major 
findings of three model experiments are presented below. 

Figure 1 summarizes the results of experiment under o clear sky. The 
model simulated PBL height, 0 at the surface and 25-m level, heat and 
momentum fluxes in the surface layer, and the eddy conductivity at the 
25-m level are shown as a function of local time for a 24-h period starting 
at 6 am (0 h on the figure). The initial PBL height was set at 375 m. During 
the first 12-h, solar heating resulted in a large amount of upward heat flux. 
The unstable PBL grew to about 800 m at the noon and reached a maximum 
of about 1050 m near 6 pm The maximum air temperature lagged behind the 
largest temperature difference between the surface and air as well as the 
maximum upward heat flux by about 6 h. The diurnal range for the surface 
temperature was close to 25°C and for the air temperature was about 4°C. 

in tne period after sunset and before sunrise the stable environment 
prevailed Air was warmer than the surface and slightly downward heat 
flux was observed. The stable PBL was much shallower than the unstable 
one There was a sharp drop in the height and the eddy conductivity during 
the transition period from the unstable to stable PBL but a relatively minor 
decrease in the momentum flux took place. This simulated the collapse of 
unstable PBL after the cut-off of the surface heat source and the slow 
re-development of stable PBL at night. 

Figure 2 is similar to Fig. 1 except for the experiment under a cloudy 
sky A layer of low cloud was placed between the 1.25 km and 2 km model 
levels Comparisons between Figs. I and 2 reveal that 

1) The model cloud decreased the unstable PBL height in the daytime, 
wmie incr eased the stable PBL height slightly at night. The vertical mixing 
in the unstable PBL was suppressed as indicated by the lower K value. 

2) The model cloud had modulated the surface temperature considerably 
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resulting in much smaller diurnal range. Also, it caused few degree cooling 
of air in late afternoon and slightly warming of air in the early morning. 

3) The model cloud reduced the temperature differences between the 
surface and air and consequently the heat fluxes, but had relatively minor 
impact on the momentum fluxes. 

Figure 3 shows the results from the third experiment. In this 
experiment, the model integration started at the local noon instead of 6 am 
and the initial PBL height was set at 750 m. The main purpose of this 
experiment is to see how the model behaves during the transition periods 
from unstable to stable then back to unstable air. Also, it will reveal how 
sensitive the model is to the initial PBL height. Comparing Figs. 3 with 1, 
we find that except for the first few hours the model behaved very similarly 
between the two experiments. For examples, in both cases the PBL grew to 
slightly over 1 km at 6 pm, and the diurnal ranges in the surface and air 
temperatures in the two experiments were very close. 


6. CONCLUSIONS AND RECOMMENDATIONS 

An one-dimensional model coupling the PBL processes and atmospheric 
radiation has been developed. The simulated results indicate that the model 
is capable of producing physically realistic solutions. A few more 
sensitivity tests with different soundings and surface parameters will be 
desizable for further understanding. For example, we may alter the 
roughness length or the Bowen ratio and see how the one-dimensional 
system responds to the changes. We hope to obtain some PBL observations 
for detailed model verifications in the near future. Comparisons with other 
approaches in modeling the PBL processes will also be very interesting. 
After carefully assessment of its performance under various conditions, the 
PBL-radiation scheme will then be incorporated into the LAMPS. 
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Table 1 


z 

(km) 

16.0 

M.O 

12.0 

10.5 

9.0 

7.5 

6.0 

1.5 

3.0 

2.0 

1 .25 

0.75 

0.375 

0.25 

0,0 


Observed data used in the numerical experiments 


u 

V 

(m/s) 

(m/s) 

10. 

-5.8 

32,2 

3.0 

15.2 

20.9 

15.9 

16.3 

21 .7 

8.0 

18.3 

12.0 

18.9 

11.1 

7.8 

-10.0 

-3.1 

-5.9 

-5.6 

-7.3 

-6.0 

-8.5 

-1.7 

-9.1 

-1.1 

-6.0 

2.5 

-3.0 

2.5 

-3.0 


6 

q 

(K) 

(g/kg) 

103 

0.01 

359 

0.05 

311 

0.1 

337 

0.2 

332 

0.3 

327 

0.3 

322 

2.0 

316 

2.8 

310 

3.3 

301 

6 

299 

9 

297 

10 

296 

11 

295 

12 

295 

12 
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Fig. 1 Model simulated (a) the PBL height, (b) heat flux, (c) ground 
(dashed line) and 25-m air temperatures, (d) eddy conductivity, and (e) 
momentum flux as a function of local time starting at 6 am for the case of 
clear sky. 
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